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Abstract 

To seek for a possible origin of fractal pattern in nature, we perform a molecular 
dynamics simulation for a fragmentation of an infinite fee lattice. The fragmenta- 
tion is induced by the initial condition of the model that the lattice particles have 
the Hubble- type radial expansion velocities. As time proceeds, the average density 
decreases and density fluctuation develops. By using the box counting method, it is 
found that the frequency-size plot of the density follows instantaneously a universal 
power-law for each Hubble constant up to the size of a cross-over. This cross-over 
size corresponds to the maximum size of fluctuation and is found to obey a dynam- 
ical scaling law as a function of time. This instantaneous generation of a nascent 
fractal is purely of dynamical origin and it shows us a new formation mechanism of 
a fractal patterns different from the traditional criticality concept. 
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Since the pioneering work of Mandelbrot [1] a large effort has been directed 
towards understanding the ubiquitous manifestation of fractal structure in 
nature. Since fractal structure is intimately related to a scaling property, a 
customary way of identifying the fractal is to check if the power-law relation 
of the frequency as functions of the size holds. Malcai et al. [2] in this way 
have listed the fractal dimensions and the corresponding scaling ranges for a 
variety of physical systems. They pointed out that the finiteness of the fractal 
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range is a typical feature for physical systems which is in sharp contrast with 
an ideal mathematical concept of fractal. 

Let us therefore settle our target as "why and how a power-law relation in 
frequency-size relation appears frequently in nature no matter what its range" . 
It seems natural to adopt some fractal growth models to challenge this subject 
or to adopt some kinds of cellular-automata approaches. In these models, the 
power-law in frequency-size relation is established at an asymptotic time. Now 
the question is if we can detect a bud of fractal pattern in the short-time regime 
long before the asymptotic time. Anatomy of this formation process at short 
time is the subject of this paper. 

In a typical fractal growth model of diffusion limited aggregation (DLA) [3], 
the short-time behavior is not studied directly. Its behavior is studied in the 
simulation study [4] and experimental study [5] of kinetic roughing model, in 
the model for fracture process [6] and in several cellular-automata approaches 
used to study the self-organized criticality [7,8]. These models in short-time 
regime are, however, concentrated on the phenomena of dynamical scaling 
relation [9] where the time-dependence of some growth parameter is inves- 
tigated. The behavior of the growth parameter, however, is not an sufficient 
answer to the growth mechanism of fractal pattern asymptotically achieved. 

In order to investigate the short-time behavior of the formation of fractal 
pattern, we propose a new simulation model for the fragmentation of an infinite 
fee lattice. For this purpose, we adopt a box-counting method [1] to search 
for the fragmentation process instead of the fragment mass distribution. It is 
because the mass distribution is only defined for the low-density limit where 
a fragment is separated each other, which is not the case for the short-time 
regime where the density is not low enough to identify the fragments. In fact, 
we will see later that the box counting method is an appropriate method to 
investigate the bud of fractal, which starts as far as the system starts to expand 
from its saturation density. 

The dynamics of the fragmentation is caused in our model by the Hubble-like 
radial expansion that is imposed as the initial condition of the molecular dy- 
namics simulation [10,11,12,13]. The Hubble expansion is the simplest and 
most symmetric driving force for the fragmentation since any point in the 
system is equivalent to any other point, a perfectly local equivalence for ev- 
erywhere when we deal with an infinite system. Infinite system is chosen to 
eliminate the effect of surface and thus we concentrate on a fragmentation 
initiated by the fracture occurring in the volume. We perform a molecular 
dynamics simulation using the Lennard- Jones (L-J) potential with a cut-off 
length r c = 2.5 [13]. As the initial condition, the particles are located at the 
lattice points of the fee lattice in the ground state, however, 10% of randomly 
selected lattice points are removed so that the initial density p becomes 0.95 in 
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a reduced unit (all quantities in what follow are measured in the reduced unit). 
This imperfection is introduced to prepare enough number of different initial 
states of the fragmentation. Our main results are essentially independent of 
the vacancy fraction, in the range 5% to 15%. The initial velocity of a particle 
at r is given by hr where h is a dimensionless parameter. Thus the system 
mimics a uniformly expanding Hubble universe. When solving the equation of 
motion for the infinite system represented by one primitive cell and 26 replica 
cells, we impose a generalized periodic boundary condition [10,11,12,13], un- 
der which a replica cell with the origin at Rj moves with the velocity foRj. This 
means that the time evolution is completely determined by the one parameter 
h. 

A complete nascent fractal structures emerge only for suitably selected value 
of h. When h is too small, the effect of the attractive force dominates the 
expansion and one single droplet whose mass equals to the total number of 
particles in a unit cell is formed. Since this structure depends on the unit cell 
size, we judge it unphysical. When h is too large, the expansion dominates the 
response of the system and the solutions to the equation of motion result in 
a simple enlarged copy of the initial configuration. Since the initial configu- 
ration with near saturation density is almost uniform, the resultant structure 
has trivial dimension 3, even after the expansion. It is only with h in an in- 
termediate range, namely 0.1 < h < 0.4 that our numerical simulations result 
in an emergence of intriguing fractal structures. 

In the box-counting method, a unit cubic cell of the size R ce \\ of expanding 
matter is divided into smaller cubic boxes of the size S = R cc n/n where n 
is a positive integer. We count the number N(d) of small boxes that contain 
at least one particle. The fractal dimension Df (box counting dimension) is 
introduced with the relation N(5) = Rl cn /5 Df , which is rewritten log[iV(5)] = 
31og[-R cc n] — Df\og[5]. When \og[N(5)] is proportional to log [5], Df becomes 
constant and the relation between N(S) and 5 satisfies a power-law. In the 
following discussion, we will use the log — log plot of N(S) and 5 to obtain Df. 

Figure 1 shows the time evolution of snapshots of the fractal structures and the 
results of the box counting method applied to these structures. The left column 
corresponds to h = 0.1 and the right column to h = 0.2. For each column, 
the time evolution proceeds from the top to the bottom. As time proceeds the 
density decreases. The three graphs presented are for the average densities 
p =0.20, 0.10, and 0.05. For our studies, we prepare 50 initial configurations 
each of which is composed of 1,235 particles (fee lattice with 10% vacancies). 
N(S) is calculated as an average of these 50 configurations. Note that for 
5 < 1, N(S) is always equal to the total number of particles because two 
particles cannot come closer than 8 ~ 1 due to the repulsive term of the L- 
J potential. From lines fitted by the method of least squares and from the 
location of the intersection point between the two lines we conclude that the 



3 



physical size of the fractal structure grows with decreasing density and that 
Df does not change significantly during the time evolution for a given h, that 
is, 1.95 < D f < 1.98 for h = 0.1 and 1.66 < D f < 1.74 for h = 0.2. The value 
of Df for h = 0.1 is larger than that for h = 0.2 and the range of scale where 
Df remains constant for h — 0.1 is wider than for h = 0.2. We checked if these 
results are independent of the number of particles in the unit cell and found 
that a cell composed of 9,878 particles leads to almost the same results. In the 
snapshot picture, we show an example of a configuration with 9,878 particles. 

A similar behavior is observed in the case of more rapid expansion. In Fig.2, 
we show snapshots of the fractal evolution and the corresponding plot of the 
results of the box counting method for h = 0.3 and h = 0.4 at p = 0.05. From 
Figs.l and 2, we conclude that a larger h gives a smaller Df and narrower range 
of scale. Note that the intersection point of the 3-dimensional and Df lines 
corresponds to the maximum size of the density fluctuation in the system. 
Beyond this crossover point, the distribution of matter is uniform and its 
dimensionality is three. 

A summary of fractal structure for h — 0.1 and h = 0.2 is shown in Fig. 3 
where the fitted lines used to determine Df in Fig.l are depicted together 
with additional lines for p = 0.40. For each figure, the four lines corresponding 
to p = 0.40,0.20,0.10,0.05 are almost overlapping for small values of 5. The 
intersection points of these lines with the four non-overlapping 3-dimensional 
lines in Fig. 3a and Fig3b are shifted toward larger 5 values with decreasing 
density. The overlap of the four lines for small value of 5 is the most important 
result of this study. It means that there is a nontrivial physical process that 
tunes the box counting so regularly. This feature is in sharp contrast with the 
results for h — > oo as shown in fig.3c. In this case, the effect of expansion 
dominates the system response and the density snapshot is a trivial enlarged 
copy of the initial configuration that simply results in a parallel shift of the 
3-dimensional lines in the figure. It is remarkable that the fractal structures 
born already at densities near p = 0.4 are maintained at lower densities. 
The approximate constancy of Df and the overlap of Df lines as functions 
of time clearly shows that there is a special correlations governing the whole 
system. Although we have not yet fully understood the mechanism behind this 
phenomenon, it is clearly caused by the critical balance between the externally 
imposed expanding motion and the internal potential. 

When we follow the time evolution beyond p = 0.05, the least square fit of 
two lines to the data becomes worse. This degradation occurs when the inter- 
nal expansion of connected aggregations of particles almost stops. After this 
time, only the space between the fragments expands. This is another reason 
that we used the word critical balance since we need an on-going competition 
between the expanding external field and internal system forces in order to 
get persistent fractal behavior. As a result, the range of scale for which we 
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get a good least-square fit for the fractal dimension Df is limited to about 
one-order of magnitude. 

Finally, we will see the meaning of the cross-over point. As was discussed 
before, the box-size of this crossover point is a measure of the maximum size 
of the density fluctuation at a given time. The time- dependence of this size is 
shown in Fig.4, where the cross-over box-size S cross for h — 0.1 and for h = 0.2 
are plotted together with the cell size i? ce ii as a function of the quantity 1 + ht. 
The calculated points are shown by various symbols corresponding to fours 
system densities, p =0.4,0.2,0.1 and 0.05 from the left to right, and lines are 
fitting to these four points. The linear dependence of R ce ii on the abscissa is 
the original assumption of boundary motion. On the other hand, the fitting 
to cross-over size yields two curves 



These power-law dependence of the quantity on time is known as the dynam- 
ical scaling phenomena [9,4,6,5] M The dynamical exponent depends sensi- 
tively on the expansion speed and is larger for slower expansion. Their nature 
seems to be quite different from the corresponding quantities examined so far 
[9,4,6,5]. We like to stress here that when we extrapolate to t=0, the value of 
Across becomes about 1.2, which nearly coincides with the smallest size of the 
system fluctuation. We can say, in other word, that the maximum size of the 
fluctuation starts to increase in coincidence with the start of expansion. The 
bud of a fractal starts to grow from the saturation density already. 

In conclusion, we have shown that a set of many-body coupled Newton equa- 
tions together with an externally imposed initial condition of uniform expan- 
sion for an fee lattice directly leads to the birth and evolution of fractal struc- 
tures for specific ranges of the initial expansion rate. The fractal dimension 
obtained with the box counting method remains constant up to a cross-over 
scale , which obeys the dynamical scaling law as a function of time. The most 
important finding is that the nascent fractal is generated from the very begin- 
ning of the expansion, which might give a hint to the original motivation of 
this study, "why the fractal-like structure prevails in nature?" . 




for h = 0.1, 
for h = 0.2. 
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Fig. 1. Snapshots [14] of fractal structures and log- log plots of results of the box 
counting method for six sets of h and p values. The left column is for h = 0.1 and 
the right for h = 0.2. For each column, the time evolution proceeds from the top 
(p = 0.2) to the bottom (p = 0.05) where p stands for the average density. In the 
log-log plot, 5 is the size of the box used for counting and N(5) is the number of 
boxes that contain at least one particle. Also plotted are vertical bars corresponding 
to the variance in the 50 sample cases. They are nearly the size of the dots. The 
lines in the log-log plot are obtained from a least square fit, the slopes of which are 
fractal dimensions. The images are drawn rescaled to the same size. The size before 
rescaling were J R cell =36.7 (p = 0.20), 46.2 (p = 0.10) and 58.2 (p = 0.05). 
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Fig. 2. Snapshots [14] of fractal structures and log- log plots of results of the box 
counting method at average density p = 0.05 for h = 0.3 (top) and h = 0.4 (bottom). 
Various quantities plotted are the same as those of the bottom of Fig.l. 





Fig. 3. Summary of fractal structure growth in expanding matter. The top two 
figures show the lines appearing in Fig.l together with additional lines for p = 0.40. 
The bottom figure is for h — > oo in which case the dimensionality of the space is 3 
at any scale, a typical feature expected for the enlarged copies. 
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Fig. 4. The size of the crossover point 6 cr0 ss together with the cell size i? cc n as a 
function of the quantity 1 + ht where t is the elapsed time. The calculated values 
are denoted by symbols corresponding to p = 0.4,0.2,0.1,0.05 from the left to the 
right. Lines are fitting to these four symbols which expression are also given in the 
figure. 
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